Modeling Uncertainty in Pith Offsets

Setup

library(matrixStats)
library(dplyr)
library(forcats)
library(rstanarm)
library(bayesplot)
library(shinystan)
library(ggplot2)
library(ggridges)
library(here)
options(mc.cores = parallel::detectCores(logical = FALSE) - 1)
if(.Platform$OS.type == "windows") options(device = windows)

Overview

The second component of ageing error arises from cores that missed the chronological center of the stem. In such cases it is necessary to estimate the pith offset, or the number of rings between the innermost observable ring and the pith. Duncan (1989) developed a geometric formula to calculate the pith offset, assuming concentric circular rings. The model first calculates the missing radius \(r\) given the measured length \(L\) and height \(h\) of the largest fully visible arc formed by the innermost growth rings:

\[ r = \dfrac{L^2}{8h} + \dfrac{h}{2}. \]

The mean width of the innermost \(N_{d}\) visible rings, \(\bar{d} = (1 / N_d) \sum_{i=1}^{N_{d}}{d_i}\), is then used as an estimate of the mean unobserved ring width to calculate missing rings as \(r / \bar{d}\).

Propagating uncertainty blahblah…

Data Collection

Describe sampling protocol, criteria for excluding cores, etc…

# 3 innermost ring widths for pithless-cored trees in all plots
inner_rings_raw <- read.csv(here("data","PithlessTCs_InnerRingWidths.csv"), 
                            header = TRUE, stringsAsFactors = FALSE)
inner_rings <- inner_rings_raw %>% 
  rename(patch = Patch, tree = Core_ID, ring = Ring_ID, width = Ring.Width) %>% 
  mutate(ring = substring(ring, nchar(ring)))

# arc length (L) and height (h) of first full ring
# remove trees not present in inner_rings, those with unknown plot ("RLARGE01", "RLARGE02")
# and those in plot "R0X" ("R0XT01a")
duncan_raw <- read.csv(here("data","Duncan_estimates.csv"), header = TRUE, 
                       na.strings = c("NA","#DIV/0!"), stringsAsFactors = FALSE)
duncan <- duncan_raw %>% rename(tree = `Core..`, h = H, width_1st_full_ring = X1st.full.ring,
                                r = Missing.Distance.To.Pith..DTP.,
                                width_inner3_rings = Width.innermost.3.rings..WI3R.,
                                mean_width_inner3_rings = Mean.WI3R,
                                rings_to_pith = Estimated...rings.to.pith..RTP.,
                                subtract_if_not_arching = If.last.vis.ring.not.arching..extra.rings.to.subtract,
                                rings_to_pith_adj = Adj.rings.to.pith,
                                duplicate = Adj.rings.to.pith.values,
                                downward_curvature = downward.curvature) %>% 
  mutate(patch = inner_rings$patch[match(tree, inner_rings$tree)], .before = tree) %>% 
  select(-duplicate) %>% filter(!is.na(patch) & !(tree %in% c("R0XT01a", "RLARGE01", "RLARGE02"))) 

Inspect data objects…

Hierarchical Linear Model

Describe modeling approach…

\[ \begin{aligned} \text{log}(d_i) &= \mu_{\text{tree}[i]} + \epsilon_i = \beta_{\text{patch}[i]} + b_{\text{tree}[i]} + \epsilon_i \\ b_{\text{tree}} &\sim N(0,\sigma_\text{tree}) \\ \epsilon_i &\sim N(0,\sigma_d). \end{aligned} \]

duncan_lmer <- stan_lmer(log(width) ~ patch + (1 | tree), data = inner_rings, 
                         chains = getOption("mc.cores"), iter = 5000, warmup = 1000)

print(duncan_lmer,2)

stan_lmer
 family:       gaussian [identity]
 formula:      log(width) ~ patch + (1 | tree)
 observations: 366
------
                   Median MAD_SD
(Intercept)        -2.01   0.12 
patchBurned Forest  0.66   0.18 
patchForest         0.57   0.16 
patchTransition     1.52   0.17 

Auxiliary parameter(s):
      Median MAD_SD
sigma 0.28   0.01  

Error terms:
 Groups   Name        Std.Dev.
 tree     (Intercept) 0.634   
 Residual             0.281   
Num. levels: tree 122 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

Posterior Predictive Checking

Discuss posterior predictive checking…

yrep <- posterior_predict(duncan_lmer)
indx <- sample(nrow(yrep), 2000)

Posterior predictive marginal density…

Posterior predictive check for mean log ring width…

Posterior predictive check for SD of log ring width…

Normal QQ plot of tree-level intercepts, grouped by patch…

Normal QQ plot of observation-level residuals, grouped by patch…

Pith Offsets

We transformed the posterior samples to calculate the posterior distribution of median inner ring width for each tree, \(\widehat{d_{\text{tree}}} = \text{exp}(\mu_{\text{tree}})\), which we then used in place of the sample mean \(\bar{d}\) to estimate the missing rings to pith as \(r / \widehat{d_{\text{tree}}}\). We used the median rather than the mean as a more robust measure of location for the lognormal data distribution, although in this case the difference is slight.

d_hat <- exp(posterior_linpred(duncan_lmer, newdata = duncan))
rtp <- sweep(1/d_hat, 2, duncan$r, "*")  # keep continuous version for plotting

# Attach posterior draws to data
rings_to_pith <- data.frame(duncan[rep(1:nrow(duncan), each = nrow(rtp)), c("patch","tree")],
                            iter = rep(1:nrow(rtp), nrow(duncan)), 
                            rings_to_pith = as.vector(rtp))
rownames(rings_to_pith) <- 1:nrow(rings_to_pith)